### Replication file for Bas & Schub "Peaceful Uncertainty"
### (I)  Generate Figures 4 and 5
### (II) RAND Benchmarking from Supporting Files


##### (I) Generate Figures 4 and 5 #####
rm(list=ls())
library(foreign)

# set working directory
dat<-read.dta("BasSchub_ISQ_Shocks.dta")

######
## Figure 4: Left Panel

# subset iraq-iran data
irir<-dat[dat$initiator==645 & dat$target==630 & dat$year>1974 & dat$year<1981,]

# plot
par(mfrow=c(1,1))
par(mar=c(3,3,2,1))
par(oma=c(1,1,0,1))
plot(irir$year[irir$year<1978],irir$Pactual[irir$year<1978],pch=19,ylim=c(0.15,0.65),ylab="",xlab="",xlim=c(1975,1980),cex.axis=1.3,cex=1.2)
#points(irir$year,irir$Pmean,type="l",lwd=3)
points(irir$year,irir$Pup,type="l",lwd=2,lty=3)
points(irir$year,irir$Pdo,type="l",lwd=2,lty=3)
points(irir$year[irir$year>1977],irir$Pactual[irir$year>1977],pch=0,col="black",cex=1.2)
points(irir$year[irir$year>1979],irir$Pactual[irir$year>1979],pch=15,col="grey60",cex=1.2)
mtext(cex=1.4,side=2,outer=FALSE,line=2.4,"Iraq Share of Dyadic Capabilities")
mtext(cex=1.5,side=3,outer=FALSE,line=0.7,"Iraq-Iran: Projections and Realizations")
text(1975.5,0.32,"cutpoints",cex=1.2)
segments(1975.5,0.31,1975.7,0.271)
text(1977,0.4,"realized value",cex=1.2)
segments(1977.4,0.39,1977.9,0.367)
legend("topleft",col=c("grey60","black","black"),pch=c(15,0,19),legend=c("shock, neutral future trend","shock, favorable future trend","no shock"),bty="n",cex=1.3)


######
## Figure 4: Right Panel

## substantive quantities calculated in "BasSchub_Shocks_ISQ_Main.do"
range<-c(1,2)
est<-c(1.79,3.10)
lo<-c(1.23,1.93)
hi<-c(2.45,4.76)

# plot
par(mfrow=c(1,1))
par(mar=c(3,3,3,1))
par(oma=c(1,1,0,1))
plot(range,est,pch=20,cex=1.6,ylim=c(-0.2,5),xlim=c(0.5,2.5),xaxt="n",xlab="",yaxt="n")
segments(range,lo,range,hi,lwd=3,col="grey50")
points(range,est,pch=20,cex=1.8)
abline(h=0,lty=2)
mtext(cex=1.6,side=3,outer=TRUE,line=-1.5,"Effect of Power Shocks")
mtext(cex=1.4,side=3,outer=TRUE,line=-3,"Iraq-Iran Dyad")
mtext(cex=1.3,side=2,outer=FALSE,line=2.5,"Predicted Probability of War (%)")
axis(1,at=c(1:2),labels=c("1977: No Shock","1980: Shock"),cex.axis=1.3,tck=0)
axis(2,at=seq(0,5,by=1),cex.axis=1.3)



######
## Figure 5: Left Panel

## substantive quantities calculated in "BasSchub_Shocks_ISQ_Main.do"

range<-seq(-0.04,0.04,by=0.01)
est<-c(0.60,0.69,0.78,0.88,0.98,1.08,1.18,1.28,1.39)
lo<-c(-0.37,-0.28,-0.19,-0.08,0.02,0.14,0.24,0.35,0.45)
hi<-c(1.80,1.78,1.87,1.94,2.04,2.15,2.27,2.35,2.47)

rivs<-dat[dat$rival==1,]

# plot
par(mar=c(3,3,3,1))
par(oma=c(1,1,0,1))
plot(range,est,pch=20,cex=1.6,ylim=c(-0.8,2.6),xlim=c(-0.045,0.045),type="l",lwd=4,xlab="",ylab="",xaxt="n",yaxt="n")
points(range,lo,range,type="l",lty=2,col="grey50",lwd=2)
points(range,hi,range,type="l",lty=2,col="grey50",lwd=2)
abline(h=0,lty=3,lwd=2)
mtext(cex=1.6,side=3,outer=TRUE,line=-2,"Conditional Effect of Power Shocks")
mtext(cex=1.3,side=2,outer=FALSE,line=2.5,"Percentage Point Increase in Pr(War)")
axis(1,at=c(-0.03,0,0.03),labels=c("Favorable","Neutral","Unfavorable"),cex.axis=1.1,tck=0)
mtext(cex=1.3,side=1,outer=FALSE,line=2.4,"Estimated Trend")
axis(2,at=seq(-0.5,2.5,by=0.5))
rug(rivs$powershift[seq(1,length(rivs$powershift),by=7)],lwd=0.2,col="grey50")


######
## Figure 5: Right Panel

rjs<-dat[dat$rival==1 & dat$Pmean>0,]
rjs$Pabsshock<-rjs$Prealizeddiff/rjs$Pmean
rjs<-rjs[complete.cases(rjs$Pabsshock),] 

cuts<-c(seq(-0.6,0.5,by=0.1),1)
mat<-matrix(NA,nrow=(length(cuts)-1),ncol=6)
for(i in 1:(length(cuts)-1)){
lb<-cuts[i]
ub<-cuts[i+1]
grab<-rjs[rjs$Pshiftactual>=lb & rjs$Pshiftactual<ub,]
warnoshock<-mean(grab$war2[grab$Pshock==0])
warshock<-mean(grab$war2[grab$Pshock==1])
mat[i,1]<-lb
mat[i,2]<-ub
mat[i,3]<-lb+0.05
mat[i,4]<-warnoshock
mat[i,5]<-warshock
mat[i,6]<-nrow(grab)
}
mats<-mat[7:12,]
df<-as.data.frame(mats)
colnames(df)<-c("lb","ub","mid","warnoshock","warshock","grab")

# plot
par(mar=c(3,3,3,1))
par(oma=c(1,1,0,1))
plot(df$mid,df$warshock,col="grey60",pch=18,xlim=c(0,0.6),ylim=c(0,0.2),xaxt="n",xlab="",yaxt="n",ylab="",cex=1.5)
points(df$mid,df$warnoshock,col="black",pch=19,cex=1.2)
rug(rjs$Pshiftactual[rjs$Pshock==1 & rjs$Pshiftactual>0],col="grey60",lwd=0.3,side=3)
rug(rjs$Pshiftactual[rjs$Pshock==0 & rjs$Pshiftactual>0],col="black",lwd=0.3)
abline(v=seq(0,0.6,by=0.1),lty=3)
mtext(cex=1.5,side=3,outer=TRUE,line=-1.8,"Year-over-Year Power Shifts:")
mtext(cex=1.4,side=3,outer=TRUE,line=-2.9,"Shocks vs. No Shocks")
axis(1,at=seq(0.0,0.6,by=0.1),cex.axis=1.2)
axis(2,at=seq(0.00,0.20,by=0.05),cex.axis=1.2)
mtext(cex=1.3,side=1,outer=FALSE,line=2.4,"Annual Shift in Power Balance")
mtext(cex=1.3,side=2,outer=FALSE,line=2.4,"Probability of War")
text(0.34,0.13,"shock",col="black",cex=1.2)
text(0.37,0.02,"no shock",col="black",cex=1.2)
segments(0.32,0.124,0.345,0.115)
segments(0.385,0.014,0.355,0.004)




##### (II) RAND Benchmarking from Supporting Files #####

######
## RAND Benchmarking and Figure A2
library(foreign)
rand<-read.dta("BasSchub_ISQ_Shocks_randbenchmarking.dta")

# plot
plot(rand$estmilex_A,rand$implied1990v2,pch=20,,yaxs="i",xaxs="i",ylim=c(-0.1,41),xlim=c(-0.1,35),ylab="RAND Projection 1990",xlab="Our Projection",main="Projected Military Spending")
segments(0,0,40,40,lwd=3)
text(9,27,"R2=0.84")

# Scaled model: R2=0.84
summary(lm(implied1990v2 ~ estmilex_A,data=rand))

# Unscaled model: R2=0.88
summary(lm(rand ~ estmilex_A,data=rand))






